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Abstract. In this paper we discuss an abstract iteration scheme for the calcu- 
lation of the smallest eigenvalue of an elliptic operator eigenvalue problem. A 
short and geometric proof based on the preconditioned inverse iteration (PIN- 
VIT) for matrices [Knyazev and Neymeyr, (2009)] is extended to the case of 
operators. We show that convergence is retained up to any tolerance if one only 
uses approximate applications of operators which leads to the perturbed precon- 
ditioned inverse iteration (PPINVIT). We then analyze the Besov regularity of 
the eigenfunctions of the Poisson eigenvalue problem on a polygonal domain, 
showing the advantage of an adaptive solver to uniform refinement when using 
a stable wavelet base. A numerical example for PPINVIT, applied to the model 
problem on the L-shaped domain, is shown to reproduce the predicted behaviour. 



1. Introduction 

In problems arising from physics and engineering one is interested in finding the 
smallest eigenvalue and/or corresponding eigenfunction of a given elliptic partial 
differential equation. Depending on the context, this can be for example the lowest 
vibrational mode in mechanics, or the ground state energy in chemical structure 
calculation. 

In a standard way an eigenvalue problem is posed in a weak formulation IH. We are 
looking for the smallest eigenvalue A € M and corresponding eigenvector u ^ V, 
such that 

(1) a{u, v) = A(m, v), for all v ^V, 

where V is an appropriate Banach space (e.g. H^) that is a dense and continuously 
embedded subspace of a Hilbert space H (e.g. L2) with inner product (•, •). We 
assume that a is a bounded, symmetric and strongly positive bilinear form. Fur- 
thermore we assume that the smallest eigenvalue Ai is simple and well separated 
from the rest of the spectrum. 

Using Finite Element Methods (FEM), the eigenvalue problem can be efficiently 
solved numerically. However when the eigenfunction exhibits singularities, one 
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has to use adaptive strategies to retain efficiency. For the iterative mesh generation 
one often uses local error estimators or indicators ||23l or, more recently, also dual 
weighted residual based goal oriented error estimators [211181. In practice these 
methods perform well, but optimal convergence rates cannot be proven yet lilSil . 

A benchmark for optimal convergence rate is the nonlinear best A^-term approxi- 
mation of the solution [12J. Therefore one can expect an adaptive algorithm at best 
to calculate an approximation to the solution with an effort which is proportional 
to the degrees of freedom needed for a best A^-term approximation of the same 
accuracy. 

In this sense, a recent article (9) following the spirit of 16] showed optimal conver- 
gence of a perturbed preconditioned inverse operation for the solution of elliptic 
eigenvalue problems . As a basis one uses the operator formulation 

(2) Au = XEu, 

where A corresponds to the bilinear form a of equation dTjl and E results from the 
//-inner product. Then given a preconditioner P for A one determines 

(3) = r;" + aP-\Av'' - fi{v'')Ev'^) 

up to an accuracy in each step, where a is an appropriate step length and fi is the 
Rayleigh quotient. The proof of convergence uses the fact that in a neighborhood 
of the eigenfunction the iteration is contracting for the part perpendicular to the 
corresponding eigenspace. 

In the case of matrices a more geometrical proof is known f20l. This proof assures 
convergence to the smallest eigenvalue for all starting vectors whose Rayleigh quo- 
tient lies between the first and the second eigenvalue. Therefore the domain of con- 
vergence can be substantially bigger compared to the alternative proof. The first 
aim is therefore to extend the proof from the matrix case fSO"] to abstract spaces. 
This will substantially shorten our proof and improve our result from |9]. More- 
over a more geometric and intuitive interpretation of the iteration is possible. We 
will show that, in order to retain convergence, each operator application has to be 
performed only with an accuracy proportional to the current error in the eigen- 
function. We will show how the idealized iteration can be performed using only 
approximate operator applications which is a common practice in adaptive wavelet 
methods 00. 

Our second aim is to apply the present abstract iteration to wavelet discretization. 
We will apply the adaptive wavelet algorithm to the model problem of a planar 
Poisson eigenvalue problem on a polygonal domain. For these ansatz spaces the 
rate of approximation of an eigenfunction is determined by the regularity of the 
function in terms of Besov spaces. Therefore we will first determine this kind of 
regularity for the eigenf unctions. It will be shown that the eigenfunctions can be 
approximated arbitrarily well provided that the wavelets have a sufficient number 
of vanishing moments. This is in contrast to Sobolev regularity, where the biggest 
inner angle of the domain restricts the smoothness of the eigenfunctions. There- 
fore, for domains with reentrant corners, the adaptive scheme is superior to uniform 
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refinement. We conclude the practical part by providing some numerical results for 
the L-shaped domain. 

We will proceed along the following line. First, in Section 2, we fix notation and 

will rewrite the eigenvalue problem in terms of operators. After that the conver- 
gence of the abstract iteration including perturbations will be shown in Section 3. 
In Section 4 we will concentrate on perturbations resulting from inexact operator 
applications. In the last section we will apply the abstract iteration to the case of 
a planar Poisson eigenvalue problem, calculate the regularity of the eigenfunctions 
and provide numerical results for the L-shaped domain. 

2. Operator formulation 

In this section we will introduce the notation, state the basic assumptions and pose 
the problem in terms of operators. This is done using the abstract setting of a 
Gelfand triple which will simpUfy the later analysis. 

For that purpose let {H, (•,•),! • |) be a separable Euclidean Hilbert space, and 
{V,\\ ■ II) a reflexive and separable Banach space such that V C H is dense and 
continuously embedded in H, i.e. 

(4) \v\ < a\\v\\ for all v gV. 

Denote by [H*, \ ■ \^) and {V*, \\ ■ ||*) the respective dual spaces of H and V. 
The dual pairing on V* and V is given by (•, •) : V* x V ^ M.. The spaces 

V C H ^ H* cV* form a Gelfand triple by identifying H* and H by the Riesz 

representation theorem. 

Assume that we are given a bilinear form a -.V xV which is bounded, sym- 
metric and strongly positive. We will consider the problem of finding the smallest 
eigenvalue and corresponding eigenvector of the weak eigenvalue problem 

a{u, v) = X{u, v) for all v eV. 

Equivalently this equation can also be written in operator form. Through the Riesz 
representation theorem, the bilinear form a uniquely determines an operator A : 

V -^V* satisfying 

a{u, v) = {Au, v) for all u, v ^ V. 

A is bounded, strongly positive, and symmetric with respect to the dual pairing 
(•, •) in the sense that 

{Av, u) = {Au, v) , for d\\v,u G V. 

Hence there exist constants (Jq and ai such that 

(5) cfqWvW'^ < {Av,v) < oi\\v\\'^ forall-ueF. 

For the formulation of the eigenvalue problem in terms of operators we introduce 
the mapping 

E:H^H*, v^{;v) 
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which is induced by the inner product (•, •) on H. For convenience we will also 
denote its restriction E\y ^ C{V,V*)hy E. 

Now an equivalent definition of a weak eigenvalue in terms of operators can be 
made. 

Definition 1. Let A : V ^ V* he. n symmetric, bounded and strongly positive 
operator. A G M is a (weak) eigenvalue if there exists a.v £ V \ {0}, such that 

(6) Av = XEv. 

Then v is called a (weak) eigenvector. The (weak) resolvent p{A) of A is given by 
all values A G M, such that Av — XEv = / is uniquely solvable for all / € H* 
and the inverse mapping is in C{H*, V). The (weak) spectrum is given by a{A) = 
M \ p{A). The Rayleigh quotient is given by 

{Ev,v} {v,v) 

We assume that the lower part of the spectrum is discrete, that is there exist eigen- 
values < Ai < . . . < Xn of possibly higher multiplicity with corresponding 
finite dimensional eigenspace 

£k = span(nfe^i,...Wfc,nJ, k = l,...,N, 

while we suppose the rest of the spectrum is bounded from below by A > X^ 
and unbounded from above. If the latter is not the case, A and E are spectrally 
equivalent and since they induce norms on V and H both spaces coincide. In this 
paper, we will restrict ourselves to the unbounded case, only noting that the other 
case can be treated with only minor modifications. 

Whenever V is compactly embedded in H, as it is the case for eigenvalue problems 
on bounded domains, the spectrum consists only of eigenvalues and the previous 
assertions are fulfilled automatically. 

The problem we will treat in the sequel can be formulated as follows: Find the 
smallest eigenvalue Ai G M and a corresponding eigenvector ui G ^\{0} such 
that 

Aui = XiEui. 



3. Perturbed Preconditioned inverse iteration for operators 



In this section we will state an iterative method for solving operator eigenvalue 
problems and show its convergence, formulated in Theorems [3] and HI For the 
construction and the analysis we can rely on methods developed for generalized 
symmetric eigenvalue problems. In particular we will use an iteration based on 
the preconditioned steepest descent of the Rayleigh quotient. These methods were 
first analyzed in 1 13, .16 ], and recent developments were achieved in I3l l2Tl l20l . 
Generalization of such iteration schemes to operators were also considered in ll24]| . 
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Let us again stress that in contrast to the above references our iteration will be 
formulated in the infinite dimensional space V, not in an associated discretized 
space. In view of a numerical realization of such an algorithm approximations are 
unavoidable in general. Therefore we will from the very beginning modify the 
Preconditioned inverse iteration (PINVIT) in allowing for a perturbation in each 
step, reflecting the finite dimensional approximation of the involved quantities. In 
Section ID we will discuss the errors resulting from this approximate application of 
operators which is common in adaptive wavelet strategies [6J. 

To state our iteration scheme, we introduce a preconditioner of the operator A, that 
is a symmetric operator P : V ^ V* , such that A and P are spectrally equivalent. 
Up to a scaling, this can be reformulated in the following way (cf. fl9l ): There 
exists a constant jp < 1 such that 

(8) \\ld-P-'A\\A<iP, 

In the case of wavelets, the discretization of P will be a diagonal matrix. 

Now we can state the basic iteration scheme. 

Definition 2. Let the starting vector £ V, ^ 0, be given and define its 
associated Rayleigh quotient as fjP = fj,{v^). A perturbed preconditioned inverse 
iteration (PPINVIT) is a sequence of vectors (v")n>o and associated Rayleigh quo- 
tients (/i")n>o generated by 

_ |~ra+l|-l~n+l 

where (C")n>o € V are perturbations. 

In order to show convergence for this sheme, we will at first generalize the results 
of II20I to the case of this perturbed operator iteration scheme: Provided that the 
perturbations are bounded by a multiple of the actual accuracy, the iteration gener- 
ates a sequence of Rayleigh quotients converging to such that the error decreases 
geometrically. Furthermore we will give a bound for the rate of convergence of the 
associated subspaces. In this context, it is obvious that the size of the perturbations 
have to match the current accuracy in the iterands to retain convergence. It turns 
out that the eigenvalue residual 

p{v) = \\Av- X{v)Ev\\a-i/\\v\\a 

is in a sense an efficient and reliable error estimator for the angle of the current 
iterand and the eigenvalue spaces. Choosing the perturbations (^")n>o in the order 
of the residuum will guarantee the convergence of the iteration. 

Theorem 3. Let v £ V, v ^ 0, such that the associated Rayleigh quotient A = 
X{v) fulfills Afc < A < Xk+i- Furthermore assume that the perturbations is 
bounded by 

U\\a/\\v\\a <JiiP(.v), where -f = jp + -f^<l. 
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Then the next step of PPINVIT ( cf. Definition |2l) with starting vector v gives v' 
and an associated Rayleigh quotient X' = A(f'), for which either X' < or 
Afc < A' < Xk+i- In the latter case 

A' — Afc 2^ \ \ X A — Afc 
T rr < q h,Xk,Xk+i)- r. 

Here q is given by 

q{l, Afc, Afc+i) = 1 - (1 - 7)(1 - Afc/Afc+i). 



Therefore, the rate of decay is only governed by the eigenvalue gap and the quality 
of the preconditioner. Note that the presence of a perturbation has the same effect 
as applying a preconditioner with a constant 7 instead of 7p. 

Besides the Rayleigh quotient one is also interested in convergence to the eigenspace, 
which is best described by the convergence of the angle between the iterand and the 
eigenspace. The following theorem states that the angle for the smallest eigenvalue 
is controlled by the magnitude of the eigenvector residual. 

Theorem 4. Let v & V, v ^ 0, such that for the associated Rayleigh quotient 
X = X{v) fulfills Ai < A < A2. Denote the angle in terms of the scalar product 
{A-, •) between v and the eigenspace by £1). Then 



■ A ( c \ ^ X{x) - Al 

'''"^^^''^^'^-^Y.-x^^y 

Moreover the eigenvector residual controls the angle, i.e. 

^-^p{v) < sm(j)A{v,Ci) < - — \ p{v) 



3A(?;)"^ ' - --^ ' - A2-A(t;)' 

Therefore convergence of the Rayleigh quotient towards Ai assures convergence of 
the angle between the iterands and the corresponding eigenspace £1 . Concerning 
the perturbation, its magnitude may be chosen proportional to the current error in 
the subspaces. 

The above statements will be proven by reducing the problem to a more simpler 
model case. 



3.1. A model case analysis. The aim of this section is to prove the convergence 
rate of a preconditioned inverse iteration for a special eigenvalue problem given by 
a bounded operator on a Hilbert space. Later on the setting of the previous part can 
be transformed to fulfill the specialized assumptions. 

Suppose that a Hilbert space X with scalar product (•, ■)x is given. Furthermore 
let i? : X — )• X be a bounded operator such that 

inf ip^=o. 

XGX\{0} {x,X)x 
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where the top of the spectrum consists of discrete eigenvalues /xjv < . . . < /ii with 
corresponding finite dimensional eigenspace X^- We are looking for the biggest 
eigenvalues /Xj of the eigenvalue problem 

Bx = fix. 

Furthermore suppose we are given an preconditioner T : X ^ X such that 

||Id — T\\x < 7t- 
Then we define the following iteration: 

(9) x' = x + ^T{Bx- fi{x)x)+T], ^(x) = -^^^^4^' 

{x,x)x 

where rj £ X is again a perturbation. The eigenvalue residual is defined by 

P{^) = l|— 7-t(-S2; - 
fi{x) 

Regarding the convergence of the sequence generated by|9]we can state the follow- 
ing estimate. 

Theorem 5. Let x X, x 0, such that for the associated Rayleigh quotient 
= ^{x) fulfills /ifc+i < A* < IJ-k- Furthermore assume that the perturbation -q is 
bounded by 

\\v\\x/\\x\\x <lr,p{x), where 7 = 7p + 7^<l. 

The above iteration step applied to a vector x then gives an output x' and an 
associated Rayleigh quotient fj,' = ^{x'), for which either p,' > p^ or fik+i ^ 
p' < /Ufc. In the latter case, 

< cr , c7=l-(l-7)- 



Proof. Neglecting the perturbation the proof of this theorem can be taken almost 
verbatim from the proof of Theorem 1 . 1 of 1201 for the matrix case, where /^min = 
0. The notation has been chosen identical to the one from 1 20 1 in order to simplify 
this transition. In the presence of perturbations, we see that the scaled iterand 
^{x)x' fulfills 

fi{x)x' = Bx — (/ — T){Bx — fi{x)x) + n{x)r]. 

Hence ^{x)x' lies in a ball with radius 7||Sx — fi{x)x\\x and center Bx since the 
distance between and Bx can be estimated by 

11(1 - T){Bx - fx{x)x) + r]\\x < ^t\\Bx - n{x)x\\x + ItjWBx - fj.{x)x\\x 

= ^\\Bx — fj.{x)x\\x- 
From then on the proof proceeds as before. □ 

The error estimation of Theorem |4] also has a counterpart in this setting. 
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Theorem 6. Let x € X, x / 0, such that for the associated Rayleigh quotient 
fx = nix) fulfills /X2 < )U < /xi. Denote the angle in terms of the inner product on 
X between x and the eigenspace Xi by (f)x{x, Xi). Then 



sin 4>x{x, ^"1) < 



/ 111 - n{x) 



fi{x) - H2 ' 

Moreover the eigenvector residual controls the angle, i.e. 

^p{x) < smcPx{x,Xi) < -/^p{x). 

Proof. Decompose x orthogonally into x = xy + xj_, where x\\ £ X\ and x_l G 
. By the definition of the Rayleigh quotient and by orthogonality 

+ Ikxilx) = l^{x)\xfx = ||x||| = llxylll + ||x_l||| 
< + M2||a;±||x. 

Then one can estimate the angle by 



sm0x = -[]-[] — < -T] — []— < \\ —r-^ ■ 

\\x\\x \\x\\\\x y M(a^) - 

For the upper bound of the second inequaUty we use the previous estimate and 
combine it with the Temple-Kato inequahty 

(/ii - /i(x))(|u(x) - H2) < \\Bx - Ai(x)x||^/||x||^, 

which will directly give the desired result. In order to prove the Temple-Kato 
inequahty let /Z = /x(x) for convenience. Then by the spectral calculus it follows 
that 

\\Bx - flx\\x = / {jj, - p,)'^d{EnX,x). 

Ja{B) 

Then the inequality follows directly from 

(yu - - Mfc) > 0, iJ,ea{B). 

For the lower estimate we see that 

p{x) < -\\Bx - Hix\\x/\\x\\x + 
We estimate the terms separately: 



\Bx - nixfx = {/J,- fii)^ d{Efj_x,x)x 

Ja(B) 

{p - pif d{Ei^x,x)x < l^i\\x-^\\x- 



f (.. ,. ^2 A(T? ™ ^ ,.2||_±||2 
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The other term gives for xi = ||x and = 1 

n{x) — Hi = {Bx, x) — {Bxi,xi) 

= {B{x — xi), (x — xi)) — /i(x)(x — xi,x — xi) 
< /xi ||x — xi Ip. 
It then follows by geometric reasoning (cf. Lemma 2 in ||9l) that 

fii- fi{x) < 2/ii||x-^||^/||x||^. 
Putting everything together gives 

p{x) < —- sm (px H ^- — sm (px < ——- sm (px, 

fi{x) n{x) fi{x) 

and hence the assertion. □ 

3.2. Reduction to the model case. The aim of this part is to show how to trans- 
form the general setting of Section |2] to the case analyzed in the previous section 

Based on the transformation Au = XEu to A~^Eu = X~^u we set 

On the space V we introduce the inner product {v, w)v = {Av, w) and the induced 
norm \\v\\v = {v, v)v = \\v\\a- Then the upper part of the spectrum of B 
consists of discrete eigenvalues /i^ = A^^ with finite dimensional eigenvalues. 
Since the spectrum of A is by assumption unbounded it follows that the infimum 
of the spectrum is zero. The Rayleigh quotient is then defined as 

{Bv,v)v {Ev,v) 

fJ.[V) = —. r = -r-. ^ 

{v^v)y {Av,V) 

As the preconditioner we set 

T = P-'^A :V ^V. 

For the quality of the preconditioner it follows due to symmetry with respect to 
that 

||Id-r||y= sup {{Id- p-'^A)v,v)v = \\ld- P~'^A\\a<-/p. 
vev\{o} 

Our aim is now to show that the iterands and the Rayleigh quotients of the iterations 
defined by PPINVIT (Definition |2l) and by the model iteration defined by Equation 
Q coincide, provided that rj = ^. For the Rayleigh quotient, it follows directly 
from equation ( 13.21 ) that = \{v)^^. For the next iterand 

v' = V ^ —T{Bv — fJ.{v)v) + T] 

H{v) 

= v + \{v)p-^ A{A-^ Ev - \{vY^v) + T] 
= v - P~^{Av - X{v)Ev) + T], 
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which is the iteration of Theorem [3] if we set r] = ^. 

For the estimation of the convergence of the Rayleigh quotients, the factor can be 
transformed to 

Afc+i — Afc _ l/nk+i — _ fJ'k — Atfc+i 
Afc+i l//^fc+i ^J'k 

For the fractions of the Rayleigh quotients 

Afc+l - A(u) - 1//^ fJ'-fJ-k+i ^i-k 

and analogously also for A'. As the factor appears on both sides the 

inequality of Theorem [3] follows from Theorem |5l 

By direct calculation it follows that 

\\-^{Bv - ^{v)v)\\v = \\Av - \{v)Ev\\a~i. 

and hence the perturbations are bounded by the same term. Moreover the upper 
and lower bounds of Theorem |4] follow immediately. 

Note that this sort of transformation has also been used in |[2ll to get rid of the 
mass matrix. 

4. Inexact operator applications 

Considering a numerical realization of the PPINVIT, the interpretation of the per- 
turbation is rather natural: In the calculation of the residual r{v) = Av — X{v)Ev, 
the corresponding coefficient vector generally has infinitely many entries, so that 
the perturbation will be related to the error of the finite dimensional approxima- 
tions of the actual residuals. In this section, we will show that the convergence of 
PPINVIT is retained as long as the approximate applications of the operators A and 
E are kept proportional to the current subspace error, which may also be measured 
by the residual through Theorem H] We will start in Section |4~T] by showing that 
inexact operator applications proportional to some e > result in an approximate 
residual r^^v) which approximates r(?;) up to a constant times e. The main diffi- 
culty here stems from the nonlinearity of the Rayleigh quotient and from the fact 
that in contrast to the finite dimensional case, we will have to deal with different 
norms for H and V. In the next section, we then devise an algorithm which de- 
termines an appropriate accuracy for the residual to preserve convergence, while 
the operator applications involved only have to be carried out proportional to the 
current subspace error. 

To start with, let us introduce the approximate operators and a measure for their 
quality. 

Definition 7. For all v £ V and e > let As{v) and Ee{v) be approximation of 
Av and Ev respectively, such that 

ll^e(v) - Av\\^, < e\\v\\, \Es{v) - £'f I* < £\v\. 
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Furthermore define the perturbed Rayleigh quotient as 



{Eeiv),v)' 



In the following we will assume that P^^ can be applied exactly, which causes 
no further restrictions: If only an approximate operator P of P is available which 
is still spectrally equivalent to A, one simply uses P instead of P in all calcu- 
lations. In the analysis, the constant jp of equation ^ is to be replaced by the 
corresponding perturbed one. In addition, the iteration may also be generalized 
by using a different preconditioner P„ in each step. The convergence results can 
be extended to such cases if the family of preconditioners {Pn)n>o is uniformly 
spectrally equivalent to A, i.e. 

\\I - P~'^A\\a < TP for all n > 0. 

In the prototype example of wavelets considered in Section [5] the discretization of 
P will lead to a diagonal matrix which can be applied exactly. 

With Definition Ul the following recursion allows a finite dimensional implemen- 
tation of the PPINVIT algorithm. 

Definition 8. Let a vector v, v 0, and a tolerance e > be given. Define the 
approximate Rayleigh quotient /i^ = ^^{v) associated to v. We let 

v' = v',-P-\A,{Ve)-fieEe{Ve)), 



Comparing this recursion with the original Definition |7] gives ^ = B^^{ri;{v) — 
r{v)), so we will have to bound the yl-norm of this expression in terms of the 
residual p{v), cf. Theorem [3] 

In the course of our analysis, it will be important to keep track of the induced error 
with respect to the given tolerance e. This is done via constants cq, ci, C2, C3 that 
will be specified in the proofs. As we are only interested in a qualitative statement, 
it suffices to know that these constants can be bounded independent of e and the 
current vector v, as long as certain requirements are fulfilled. General bounds for 
cO) ci, C2, C3 will involve constants such as a, Ai, \2,So, o"o, o"i etc. 



4.1. Inexact application of operators. First we will investigate the difference 
between the approximate and the exact Rayleigh quotient. 

Lemma 9. Let cq = min{l/2, ao}, < e < cq, v e V, v ^ 0, and Ai < fi{v) < 
A2. Then the approximate Rayleigh-quotient is bounded by 

where c\ can bounded independently of v and e. 
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Proof. Writing out the approximate Rayleigh quotient gives 

.10^ , {Av,v) + {Ae{v)-Av,v) 

^ ^'^ '~{Ev,v) + {E,{v)-Ev,v)- 

From the definition of the approximate operators and it follows that 

\{A^{v) - Av,v)\ <e\\vf < —\\v\\\, \{Ee{v) - Ev,v)\ <e\vf, 

where we used the norm equivalence of A and || • ||, equation Inserting these 
estimates in equation (flOl ) gives 

by noting that e < cq. To estimate the difference in the Rayleigh quotients, we 
subtract fi{v) to obtain 

1 + e 

Now setting 



-H{v)e < He{v) - n{v) < — i^/i(t;)e. 



ci = max I — — — —^{v), — — — —fi{v) I = — — — —fi{v) 



1 + e ' 1-e 'J 1-e 
gives |/ie(u) — fJ-{v)\ < c\e. Furthermore c\ can be bounded from above by 

ci <2(l + ao-i)A2 

since \i{v) < A2 and e < 1/2. □ 

Next we will estimate the difference between the approximated and exact residual 
with respect to a norm which will come in handy later on. 

Lemma 10. Let < e < cq, v £ V, v and Ai < fi{v) < A2. Then there exists 
a constant C2 such that for 

r{v) = Av- iJ,{v)Ev, r^iv) := A^{v) - fis{v)E^{v), 

the difference of the exact and the approximate residual can be bounded by 

\\r^{v) - r(f)|U-i/||v|U < C2e. 

Furthermore C2 can be bounded independently ofv and e. 

Proof. The norm of the difference in the residual can be estimated by 

\\re{v) - r{v)\U < \\Ae{v) - Av\\^ + ^x{v)\\Ee{v) - Ev\\^ 

+ \^J.£{v) - + \iieiv) - ii{v)\\\Ee{v) - Ev\\^, 

where we used the triangle inequality. We have to estimate H-Ei'll* and \\E£{v) — 
Ev\\^ in the norm stemming from H*. For that purpose, for / € H*, the dual 
norm can be estimated by < Using this, the definition of the 
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approximate operators (definition |7]) and the result on the approximate Rayleigh 
quotients (lemma |9l) gives 

||re(u) — r(ti) 11^, < e||f II + ^(t;)a^e||?;|| + ci£a^||t;|| + cie^Q^||u|| 

= {l + a'^[fi{v) + ci{l + e)]}\\v\\e. 

Since A^^ is bounded as a mapping between V* and V with norm cTq setting 

C2 = f^^^/^l + a'^Hv) + ci(l + e)]}, 

we can estimate 

||re(7;) - r(i;)||^-i < C2e||i'||. 

Again estimating fi{v) < A2, e < 1/2 and ci by the upper bound in lemma|9j C2 
can be bound from above independent of s and v. □ 

4.2. Error estimation from approximate residuals. In view of the convergence 
resuhs for the perturbed PINVIT as stated in Theorem |3l we may still guarantee 
convergence of the perturbed algorithm if we admit for perturbations for which 
7 = 7p + 75 < 1, where 7p < 1 is the constant entering via the preconditioner. 
For simplicity, we fix 7^ := in the sequel. To retain convergence, we have to 
bound the perturbation by the accuracy criterion 

(11) \\P~Hreiv) - r{v))\\A/\\v\\A < liP{v); 

using Lemma[TOl it is obvious that this can be guaranteed if only e is chosen small 
enough. However, approximating the residual with more accuracy than necessary 
at the present stage may lead to unnecessary costs. Therefore, this section is ded- 
icated to the analysis of an algorithm iteratively determining e{v) (for a given 
iterand v) in a way that for each step, e{v) > max(r, /9(t>)), where r is a target 
accuracy for the residual. Note that by this, a target accuracy for the subspace error 
may be fixed, cf. Theorem ID 

As a first step we define an efficient and reliable error estimator for p{v) provided 
that the tolerance e used for the computation is small enough. 

Lemma 11. Let the estimator of the residual p[v) be defined by 

Pe{v) = lke(^')l|p-i/ll^^l|p- 
For sufficiently small e it is efficient as well as reliable in the sense that 

(1 + -fp)~^Pe{v) - C2e < p{v) < (1 - -ipY'^Peiv) + 026. 

Furthermore one may choose e > p{y) such that Pe{v) ~ piv). 

Proof. The given inequalities are a direct result of the norm equivalence between 
A and P as well as Lemma [TOl 

Now choosing e < {2c2)^^ p{v) results in 

^(1 + 7Pr^Pe{v) < piv) < 2(1 - ^p)-^p,iv), 

which shows the second assertion. □ 
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Now we can use this estimator to test if the accuracy criterion ([TT]) is already met. 

Lemma 12. The preconditioned approximate residual P^'^r^{v) fulfills the accu- 
racy criterion of Equation ([TTI ) if 

(12) £ < C3Ps{v), 

where C3 is a constant independent of v and e. Furthermore e can be chosen such 
that e > p{v)- 

Proof. The idea is to bound both sided of Equation ([TTI ) and to require that the 
bounds satisfy the inequahty. Applying Lemma [TO] gives for the right hand side 

\\p-\re{v) - r{v)\\A < \\P-^A\\A\\A-\r,{v)-r{v))\\A < (1 + 7p)'/'c2e, 

while for the left hand side we simply apply the lower bound of Lemma [TT] If we 
choose 

C3 = [(1 + 7P)'/'C2 + 75C2]-'-^^ 

1 - 7P 

the assertion follows readily. 

Next we have to prove that solutions of this inequality do not get too small. In view 
of Lemma [To] choosing e < {2c2)~^ p{v) will lead to 

Pe > ^ P{v)- 

Therefore for all 

e < min((2c2)-\ C3(l - lp)/2)p{v) 
the inequality holds, showing that one can choose e > p{v). □ 

4.3. Convergence of PPINVIT. Based on the estimators of the previous subsec- 
tion a simple algorithm can be devised to calculate a suitable tolerance e: Starting 
with an initial guess we successively halve e until the accuracy criterion of Equa- 
tion ([T2I ) is met. For this e we can still guarantee that it is proportional to the actual 
error indicated by p{v). 

The last ingredient for our PPINVIT algorithm is a stopping criterion indicating 
when the exact residual has dropped below the given target accuracy r. Note that 
in this context the work load in might still become too high if the algorithm tries 
to determine according to Equation ([T2l) . while r is already smaller that r. To 
prevent this situation, we add a stopping test based on Lemma [TTI Combining this 
with the error estimator and one step of PPINVIT leads to the following algorithm: 



PPINVIT_STEP(t;,T) ^ v' 



£ := Co 
loop 

r ■= A,{v) - Pe{v)Ee{v) 

p ■= lkllp-i/ll^l|p 

if (1 — jp)~^p + C2£ < T then 
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{target accuracy reached, criterion Lemma\TTj 
return v 
end if 

if e < C3/9 then 

{accuracy for residual reached, criterion LemmaUl]! 
return v' := v — P~^r 
end if 
e := e/2 
end loop 



A similar algorithm has already been used for the determination of refined sets 
in the context of adaptive treatment of PDEs, cf. the GROW procedure in [14]. 
Note also that the intermediate computations of r^{v) need not cause much addi- 
tional computational effort; indeed, when using the APPLY algorithm from [9], 
the computation of the residual with a lower accuracy is part of the computation of 
the residual with a higher accuracy, therefore, intermediate values can be used to 
estimate the size of r{v) and to effectively approximate a suitable error tolerance 
e. 

The properties of the algorithm PPINVIT_STEP, which follow from combining 
the general theory of convergence. Theorem [3l with the results of this section are 
compiled in 

Theorem 13. For all starting vectors v ^ such that the Rayleigh quotient fulfills 
\{v) < A2, iterating PPINVIT_STEP generates a sequence of vectors, for which 
the error in the Rayleigh quotients decreases geometrically according to Theorem 
\3\where, with the above choice of^^, 7 = (1 + 7p)/2. In each step, the accuracy 
e is proportional to p{v). 

Combining Theorems \3\ and |4] shows that along with the Rayleigh quotients, the 
residual \\Av — \{v)Ev\\ji^-i will decrease. Hence the algorithm terminates after 
finitely many steps. 

For the final iterand v", there holds p{v") < t, so that in turn, sin (pAi^", £1) ^ t. 
The maximal accuracy e needed for the approximate applications of A and E stays 
bounded by 

e > max(r,p(w)). 

Combining the algorithm with a time to time coarsening worked out in more detail 
in 191 leads to an optimally convergent algorithm. It is likely that as in the case 
of boundary value problems |14| a coarsening of the updates P'^r^^v) instead of 
the iterands gives a convergent algorithm with improved performance however we 
were not able to prove this conjecture yet. 

5. Application to planar eigenvalue problems 

In this section the application of the abstract eigenvalue solver is discussed for the 
model case of a planar eigenvalue problem discretized by a stable wavelet base. 
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The corresponding adaptive procedure is superior to uniform refinement if the cor- 
responding eigenf unction has higher regularity in terms of certain Besov spaces 
than in the scale of Sobolev spaces. 

Hence first we will address the Besov regularity of the eigenfunctions and show 
that the eigenfunction can be approximated with an arbitrary high rate, provided 
that the wavelet ansatz function have sufficiently many vanishing moments. 

After that we provide some numerical tests for the case of an L-shaped domain 
with piecewise Unear wavelets. The numerical results approve the theoretical pre- 
dictions. 

As a model we introduce the Poisson eigenvalue problem on a polygonal domain 
with Dirichlet boundary condition: Let C be a bounded open polygonal 
domain with vertices x^^') , . . . , x^'^^ such that the interior angles aj at x^'^^ satisfy 
< ai < 27r for all i = 1,. . . ,d. Let the Poisson eigenvalue problem with 
homogenous Dirichlet boundary condition 

(13) —Au = Xu, onQ,, 

(14) u = 0, on dn 
be given. 

5.1. Regularity. In what follows we want to determine the regularity of the eigen- 
functions u in terms of Besov norms. To achieve this we employ the regularity 
results with respect to certain weighted Sobolev spaces as were described in [22], 
see also the references therein. 

For the construction of the weighted Sobolev spaces define for each vertex x^*) an 
infinitely differentiable cut-off function that is equal to one in a sufficiently small 
neighborhood of x*^*) and zero outside, such that the support of the functions do 
not intersect. Define := 1 — Yli=i 

Definition 14. Let the domain i7 be given as above. For the integer / > and 
/3 G M define the weighted Sobolev space ^2/3 the closure of C°°(Q) with 
respect to the norm 

j=l \a\<l 

where pi is the Euclidean distance to x^^^ and VFj(r2) are classical Sobolev spaces. 
Moreover define the weighted Sobolev space yj/jl^) closure of C^{Q) 

with respect to | K 1 1 y^; ^ (q) ■ 

Note that this definition is a slightly simplified version of the definition in para- 
graph 6.2.1 of [22]. Here we restrict ourselves to the case of polygonal domains 
and scalar /3. 

From Theorem 6.6.1 in [22,1 we can deduce the following regularity result: 



PERTURBED PRECONDITIONED INVERSE ITERATION FOR OPERATOR EVP 17 

Theorem 15. The operator — A is an isomorphism between ;_i (^) <^nd (fi) 
for all I > 2. 

Proof. In particular this is a simplified version of theorem 6.6.1 in [22]. Just note 
that f5 = I — 1 fulfills the condition 

— vr/aj < / — 1 — /3 < vr/aj, for alH = 1, . . . , d. 

since < aj < 27r for alH = 1, . . . , d. □ 

Having established the regularity result for boundary value problems with respect 
to the weighted Sobolev spaces we can use a bootstrapping technique to deduce 
the regularity of the eigenfunctions. 

Theorem 16. Let u G Hq{^) be an eigenfunction for the Poisson eigenvalue 
problem of equation (1131 ). Then 

neF2,/-i' foralH>0. 



Proof. The theorem is proven by induction. Given the existence of u as an eigen- 
function we basically use the fact that u is the solution of the boundary value prob- 
lem with right hand side \u. Hence with every application of Theorem [15] we gain 
smoothness for the eigenfunction u. 

We start the induction by noting that the eigenfunction u E L^(il). From Defini- 
tion [T4]it follows directly that 

\W\\vl^ ^ IMlL2(n)- 

Applying Theorem [15] assures that the solution u E t^i- 

Using the above fact as a starting point for an induction we will show that u G 
^2 1-1 all / = 2, 4, . . .. Suppose that u € z-i ^ even. Then by Definition 
[14] of the weighted Sobolev spaces it is obvious that u G V-ji+i- Using Theorem 



lit follows that 



' "^2,1+1 " "'^2,;+i 



and therefore u G V^2T+i- 

Hence u G V2 i-i alU = 2, 4, . . .. Noting that 

ll^ilk/'-i < \\u\\\/l 
^2,1-2 "^2,1-1 

by Definition [T4]finishes the proof. □ 



In conclusion we have established arbitrarily high Sobolev regularity for the eigen- 
functions with respect to the appropriate weights. It remains to show that this also 
implies regularity in the sense of Besov spaces. In order to show that, we use the 
norm equivalence between certain Besov norms and the weighted discrete norms 
of the coefficients of a corresponding wavelet expansion. We will follow the line of 



18 



T. ROHWEDDER, R. SCHNEIDER, AND A. ZEISER 



proof that was also used to prove Besov regularity in Q [H for the corresponding 
boundary value problem. 

For that purpose we use two-dimensional wavelets constructed from univariate or- 
thonormal Daubechies wavelets fTTI . 

Then for the set of supports 

/ = 2-^'A; + 2-^'[o, l]^ fcez^jez, 

the functions 

7]i ■=Vj,k-=^^i]{^^ --k), r/G^-, 

form an orthonormal basis in L^(M^). Here is a set of three wavelets, since the 
space dimension is two. 

For sufficiently regular wavelets there is a norm equivalence between the Besov 
norm and the discrete norm of the wavelet expansion: A function / G B^{L'^(M?)) 
for 1/r = a/2 + 1 /2 if and only if 

ii^o(/)iil2(m2) + E K/'^^)r) < 

\ve-i' iev+ J 

where P+ denotes the set of all dyadic cubes of measure < 1 and Pq is a projector 
onto a suitable subspace of L^(M^). 

From classical regularity analysis of Grisvard fTT^ we know that the eigenfunction 
u € H^/'^{Q). In order to apply the above norm equivalence, we first extend u 
to the whole space by a Whitney extension. We will denote the corresponding 
extended function also by u; it also has Sobolev regularity 3/2, i.e. u G H'^^'^(M?). 

Now we proceed along the same lines as in Q. The approximation in the interior 
of Q, and on the coarsest scale Pq (n) do not restrict the Besov regularity. All what 
remains is to estimate the wavelet coefficients in the vicinity of the vertices. For 
that purpose we introduce the distance from a fixed vertex x*-*^ as 

6j := inf llx — x*-*^ IU2, 
xeQ(i) 

where Q{I) is a cube completely containing the support of 777. Using the approxi- 
mation order of wavelets, see e.g. [5 1 we get 

\{U,VI) \ ^ 2-"^|u|,4/n(7^2(Q(7))). 

If 61 > we can further estimate for each a such that |a| = n 

Jq{i) Jq{i) I 



2,n-l ' 



hence ||n||p,/n(i2(Q(7))) < .5] 
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Next we sum up the contributions level by level. We introduce the set of indices 
corresponding to level j by 

A,:= {(/,r?),|/|=2"2i} 

and for each level appropriate layers 

A,-fc := {(/,7?) e A,-, k2-^ <6i<{k + 1)2-^}. 

The wavelets in the vicinity of the vertices x^*^ have to be treated separately and 
we restrict ourselves to the set Aj ^ for k > ki. Then it follows that 

oo oo 

k=ki {I,ri)GAj,k k=kx {I,-n)£Aj^k 

oo 

< 2-^-- E 

k=ki 

since the cardinality of Aj ^ is proportional to k. Choosing n large enough ensures 
that the sum involving k is finite. The summation over the refinement levels j then 
amounts to sum up a geometric series. 

Now the coefficients in the vicinity of the vertices can be estimated as in [8 |. Then 
finally we arrive at a regularity result in terms of Besov spaces. 

Theorem 17. Let u G Hq{^) be an eigenfunction for the Poisson eigenvalue 
problem of equation (1131 ). Then 

u G B'^{U{Vt)), for all a > 0, where 1/r = a/2 + 1/2. 

5.2. Numerical example. The last paragraph showed that the eigenfunctions of 
our model problem have arbitrarily high regularity in the sense of Besov spaces. 
However from classical regularity theory, the regularity in terms of Sobolev spaces 
is restricted by the biggest inner angle. In particular for the L-shaped domain the 
lowest eigenfunction can only be shown to be in H'^ for s < 5/3. This means 
that even for piecewise linear hat functions the convergence rate of for uniform 
refinement will be less than A^~^/^ in the norm, where N is the number of 
degrees of freedom. 

In contrast to this, usage of a piecewise linear wavelet bases with two vanishing 
moments in the adaptive algorithm PPINVIT will lead to an optimal convergence 
rate of the eigenfunction in the H^-novm with complexity N~^l'^ ||9l. Therefore, 
adaptive solution of the Poisson eigenvalue problem is superior to uniform grid 
refinement. To conclude this paper, we will demonstrate this for the test example 
of the Poisson eigenvalue problem on the L-shaped domain il = (— 1,1)^\[0, 1]^. 

The implementation of our eigenvalue solver is based on the adaptive wavelet code 
described in |25|. The wavelet basis is constructed along the lines of flOl which 
admits a diagonal preconditioner. Moreover the set of active basis functions is 
limited to the case where the successor of each wavelet is also contained in the set. 
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N (dof) 

Figure 1 . Convergence of the Rayleigh quotient for the adaptive 
algorithm with respect to the degrees of freedom. 

This resuhs in the exact evaluation of an operator application for a given index set; 
hence the Rayleigh quotient can be calculated exactly. 

As an algorithm we used one PPINVIT_STEP followed by a Galerkin eigenvalue 
solution on the fixed index set with an appropriate accuracy followed by a coars- 
ening of the iterands. It turns out that the constants from Sections 3 and 4 are too 
pessimistic and one can use much smaller values. 

In figure [T] the error in the Rayleigh quotient in the smallest eigenvalue is shown. 
As a reference value we used Ai = 9.639723844, see [4J. It can be seen that the 
error decreases like N^'^, which is as expected twice as high as the rate for the 
corresponding eigenfunction. 

Also of interest is the structure of the chosen wavelets which is shown in figure [2l 
The plot shows the center of the active ansatz functions during the sixth step for 
two different zoom levels. There one can see the self similarity in the two scales. 
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